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Abstract 

The calcium transport in biological systems is modelled as a reaction-difltusion process. Non- 
linear calcium waves are then simulated using a stochastic cellular automaton whose rules are 
derived from the corresponding coupled partial differential equations. Numerical simulations 
show self-organized criticality in the complex calcium waves and patterns. Both the stochastic 
cellular automaton approach and the equation-based simulations can predict the characteristics 
of calcium waves and complex pattern formation. The implication of locality of calcium distri- 
bution with positional information in biological systems is also discussed. 
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1 Introduction 

Many processes in living organisms such as muscle mechanics, cardiac electrophysiology, adaptation 
in photo-receptors, and cytoplasm functions involve the calcium ion transport and its physiologi- 
cal functions[l-6]. However, the exact function of Ca^+ oscillations and transport is only partially 
understood although it is believed that they involve in the intracellular communications and syn- 
chronization in the response to a local stimulus [13-15]. Even in the simplest one-dimensional case, 
the proper modelling requires many simplifications and assumptions. The nonlinearity and cross- 
coupling in the transport mechanisms usually lead to intractable governing equations. Even with 
certain approximations and simplicity, the mathematical models still lead to nonlinear reaction- 
diffusion equations if the essence of the calcium ion activities is included. 

Nonlinear reaction-diffusion systems can exhibit complex pattern formation [9-12]. The nonlinear 
system can be simulated using finite difference or finite element methods, or the alternative cellular 
automata [16]. However, existing implicit numerical solution schemes are not always robust under 
general boundary conditions. Most of the earlier work have mainly concerned the one-dimensional 
case with piecewise linear models [10]. Meanwhile, the cellular automata [17,18] have successfully 
modelled reaction-diffusion systems [16] with relative stable pattern formation. This provides a 
possibility of dealing with the difficult nonlinear problem of calcium oscillation and waves using 
finite-state cellular automata with rules derived from related partial differential equations. 

On the other hand, self-organized criticality has been found in many systems in nature [7,8] since 
its discovery by Bak [7] and his colleagues. Since calcium oscillation and waves a very complicated 
phenomenon with the characteristics of nonlinear reaction-diffusion systems, it can be expected that 
regular patterns under certain conditions. One natural question related to this is: Do the self- 
organized criticality exist in the complexity of the calcium transport? However, this is no existing 
literature addressing this question in the context of calcium transport. This is partly because the 



research on calcium activities and their mathematical modelling for biological and physiological 
processes is still at a very earlier stage [10,14-15]. 

This paper therefore first extends the existing mathematical models for nonlinear Ca^+ waves. 
Then, a new stochastic finite-state cellular automaton is then formulated to simulate the general 
nonlinear calcium transport process and pattern formation. Based on the reaction-diffusion equa- 
tions, the stochastic CA will be formulated. The pattern formation of calcium ions with realistic 
parameters will be studied. The self-organized criticality will be tested in the complex patterns of 
calcium concentration. The positional pattern formation and its implication will briefly be discussed. 



2 Mathematical Model 

2.1 Governing Equations for Calcium Transport 

In many physiological processes within cells, Ca^"*" plays an essential role in controlling cellular be- 
havior and functioning in the sense that calcium ions act as a signalling agent for a wide range of 
cellular activities. Calcium signaling is mediated through oscillation in intracellular Ca^"*" concentra- 
tion. Calcium ions can bind to a vast number of proteins and enzymes, and the binding can initiate 
a series of reactions that ends in the formation of a chemical called inositol 1,4,5-trisphosphate (IP3). 
The diffusion of Ca^"*" and IP3 through the cell cytosol can induce further release of calcium ions 
from stores in the endoplasmic reticulum (ER) through IPs-sensitive channels. These channels are 
sensitive to calcium itself, with fast activation for lower concentrations and comparatively slower 
inhibition, thus leading to the calcium-induced calcium release (CICR). Complex wave characteris- 
tics such as plane waves and spiral waves have been observed in experimental studies using Xenopus 
oocytes as pointed out by McKenzie and Sneyd [19]. However, the detail mechanisms underlying 
these oscillation and waves are only partially understood. 

There are several mathematical models for Ca^+ wave propagation in the literature, and these 
include the important models such as DeYoung and Keizer [2], Goldbeter et al. [20], Sneyd et al. [15] 
and McKenzie and Sneyd [19]. However, different cell types can results in different mathematical 
models, although the fundamentals are quite similar. The nonlinearity in the mathematical models 
can lead to complex pattern formations and spiral waves [21-24]. The model we will use in this paper 
is an extended version of a two-pool process. In the two-pool model, the calcium Ca^"*" concentration 
in the cytoplasm and the concentration in the Ca^^-sensitive pool satisfy the two-pool model [1]. 
Although it can reproduce the oscillations and waves observed in Xenopus oocytes and generating 
spiral waves in the higher dimensional situations. The parameters for our simulations are based on 
the models developed by Atri et al [1] and McKenzie and Sneyd [19]. The functioning of IP3 in the 
oscillation and waves has been reviewed in details by Sneyd et al. [25] , and th calcium signaling has 
been reviewed by Falcke [26]. 

We use the variable u{x,y,t) =[Ca^+] for the concentration (iiM) of calcium, v{x,y,t) for the 
fraction of the IP3 receptors that are active. The variable p = [IP3] can be treated as a bifurcation 
parameter for the reason discussed later. Then, the nonlinear model equations for calcium transport 
[10,13-15] can be written as 

— =DuV^u + Jf-Jp + Ji, (1) 

X^=D^W\ + g{u,v), (2) 

where J/ models the flux of calcium through the IP 3 receptor. Jp models the amount the Ca^+ 
being pumped out of the cytoplasm back into the endoplasmic reticulum or out through the plasma 
membrane. J; models the calcium leaking into the cell. We have 

Ji = a, q(u,v) = V, (4) 



and 

where k\,k2,k^,ku are constants and a,/3, 5 are parameters. In addition, D„ = is used in most 
existing models since most models base on the assumption that Ca^+ instantaneously activates the 
IP3 receptor. 

As the number of IP3 receptors remains approximately constant, thus we may assume that p is 
fixed and subsequently it can be considered as a bifurcation parameter. In fact, some studies using 
p G [0.24,0.2434] obtained very interesting results [22]. The function £(p) is the fraction of /P3 
receptors that have bound /P3 and increases as p increases. Then, the complete nonlinear model 
equations become 

X— = V, (7) 

dt kf + ^ ' 

where w is a time-dependent inactivation variable, and n — 3, d = m ^ 2. Some typical values of the 
related parameters are k = 3/iMs~^ (M=Mol/L is the molar concentration), ki — k^ — 0.7/iM, fc2 = 
0.01^M,fc^ = lfiM,ku = 0.27nM,d = 0.11, a = 0.15^lMs-'^, X = 0.2s, D„ = 20^im^s'^,-f = 
2fiMs-^ [10]. 

By using the notation s{u, v) — Jf — Jp + Ji and g{u, v) — fc™ / (/c™ + it™) — -y , we can obtain the 
steady state solution {uo,vn)- From equations (1) and (2) together with equations (6) and (7), it is 
straightforward to check that they satisfy the Turing instability conditions 

Q = ( !" !" ), (8) 

yu yv 

which requires 

tr(g) = s„ + .g„ < 0, Det(Q) = Sugv - Ugu > 0. (9) 

The model equations for intracellular calcium waves can be generally written in a system of 
reaction-diffusion equations in the form 

</)t = i?V20 + /(,/.), <^=[u vf, (10) 

where D = diag(L'u,0). The rate /(</)) is a general nonlinear function coupling the different compo- 
nents u and V. This nonlinear reaction-diffusion system can be solved using conventional numerical 
method or stochastic cellular automata. 



3 Stochastic Cellular Automaton 

Conventionally, reaction-diffusion systems can be solved numerically using finite difference or finite 
element methods. Alternatively, the nonlinear systems can be simulated by using cellular automata 
[16]. The cellular automata for reaction-diffusion systems can be formulated to correspond to the 
solutions of the related partial differential equations in a qualitative manner [4] or quantitative 
manner [5]. The macroscopic approach via cellular automata as demonstrated by [16] is always more 
efficient than explicit numerical method and can be more efficient than better numerical technique 
in many cases. 

The discretization of equation (|10p can be written as 

By taking At — Aa; — Ay ~ 1, this becomes 

C = ^M+i., + + <t>l,+i + 4>l,-i) - ^4>U + f{4>l,). (12) 
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Figure 1: Pattern formation of calcium concentrations (spiral waves on the left and rings on the 
right) using values given in section 2.1 



which can be written as a generic form 

r 

€,V= E ^>',i€+kj+i + fi<l>h), (13) 

fc,/=— r 

where the summation is over the 2r + 1 neighborhood. The coefficients Ck are determined from the 
discretization of the governing equations, and for this special case, c_i^o = c+i,o = co,-i = co,+i = 
-D, Co,o = — 4L). In fact, this is equivalent to a finite-state cellular automaton with a transition rule 
G = [gij] = 1, 2, N) from one state = [0^] = 1, 2, N) at time level n to a new 
state $"+^ = [(pij^^] = 1, 2, N) at time level n+1. That is, 

g^^. : ^» ^ 0^+1. (14) 

The basic assumption for the rule inference for stochastic automata is that the function 

5ii(Ci) = V<r(C,-,/,r,iV), (15) 

where r is a function with a range of [0,1], and V is a random variable. At every time step, a 
random number V is generated for each automaton The new state will only be updated if the 

generated random number is greater than F, otherwise, it remains unchanged. Following a similar 
derivation for ordinary differential equation [9] , the probability function F can be written as 

T{(l>2j,f, r, N) = T{e-f) =a + fee"/, (16) 

where a and b are coefficients depends on the number of finite states and other factors such as 
the accuracy of the approximation to the partial differential equations. The parameters of the 
continuum reaction-diffusion model shall be calibrated to fit the results given by the stochastic 
cellular automaton model using a least-squared procedure so as to get the related accurate transition 
rules. 



4 Simulations and Results 

By using the approach of stochastic cellular automata with finite number of states, we can now sim- 
ulate the reaction-diffusion calcium transport and nonlinear calcium waves. Numerical simulations 
are carried out on a,n N x N lattice in 2-D setting, and usually, N > 40, or up to 1500. The number 



of states is taken to be 256 in the present simulations. Different simulations with different lattice 
size are compared to ensure the simulated results are independent of the lattice size and time steps. 
In the rest of the paper, we present some results related to calcium pattern formation, nonlinear 
calcium waves and self-organized criticality of the complexity of reaction-diffusion systems. 

4.1 Pattern Formation and Calcium Waves 

From the initial random configuration, nonlinear reaction-diffusions can lead to complex patterns. 
Figure 1 shows two examples of the calcium concentrations evolving to the interesting patterns 
in 2-D configuration with 200 x 200 cells. The parameters used in the calculations are k = 
3/xMs~i (M=Mol/L is the molar concentration), ki = ks = 0.7iiM,k2 = O.Ol^M, fc„ = 1/zM, A:„ = 
0.27 fiM, S = 0.11, a = O.lbfiMs^^ , A = 0.2s. Spiral waves are formed for values D„ = 5/xm^s~^ (Fig 
la) while rings and ribbons are formed for higher value Du = 25fim'^s~^ (Fig lb). These patterns 
gradually evolve with time; however, the general characteristics of patterns only change slowly with 
time. 

The simulations imply that patterns and structures formed by local transition rules are relatively 
stable. Our present results are consistent with the other well-known studies in pattern formation 
such as [11,12] for plants and sea shells by using nonlinear reaction-diffusion equations. In addi- 
tion, our present simulations suggest that patterns arise naturally from the local interactions cither 
through rule-based/agent-based evolution in terms of relationships among the neighbourhood indi- 
viduals in stochastic cellular automata or partial differential equations in terms of system variables 
such as calcium concentrations. The initial configuration is not important. The rules of interac- 
tions or relationships between entities/individuals in the nearest neighbour are crucial factors that 
responsible for the behaviour of the system and pattern formations. 

In addition, the spatial pattern formation of calcium concentration can provide some positional 
information of calcium distribution resulting from the reaction-diffusion transport among cells. This 
may have some important implication in the mechanism of calcium functionality in biological sys- 
tems. Although many factors such as the function of proteins, genetic information and enzymes 
affect the intracellular transport of calcium, however, the calcium reaction-diffusion process itself 
could greatly contribute to the spatial distribution of calcium and thus be responsible to some extent 
for its positional information and signalling in biological systems. If it is the case, the modelling of 
calcium transport can be beneficial to the understanding of formation of the spatial structure and 
positional signalling coupled with the genetic and functional information in biological systems and 
physiological mechanisms. 

The reaction-diffusion systems of calcium transport are complex. Nonlinear waves can arise 
under certain initial and boundary conditions. For a single source, the calcium concentration starts 
to expand and form nonlinear waves as shown in the figure on the left of Figure 2 which is a 
snapshot at f = 500. Numerical simulations imply that wave speed decreases with time as expect for 
any passive diffusion systems. In addition, the for a random configuration with periodic boundary 
conditions, complex nonlinear wave pattern have observed in numerical simulations as demonstrated 
in the figure on the right (Fig. 2b) where calcium concentration varies with space and time. 

4.2 Complexity and Self-Organized Criticality 

For a lattice N = 100 x 100 with 256 states, the complexity of the cellular automata can be measured 
by its entropy S is defined as S = — X^jPilogPi, where Pi is the probability of states i [1]. For a 
finite state automaton, pi can be approximated by the calcium concentration so that 



The variation of complexity or entropy of the stochastic cellular automata with time is shown in 

Figure 3a. It is clearly seen that the complexity varies significantly at the early stage of the pattern 
formation process, then it gradually relaxes to the equilibrium at long time, indicating that the 
reaction-diffusion system is in a quasi-steady state. 




(17) 




Figure 2: Nonlinear calcium waves for a single source (left) and multiple sources (right). 
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Figure 3: Complexity and entropy variations at different time steps (left). Self-organized criticality 
in calcium pattern formation leads to the exponent 7 = 1.26 ± 0.02 (right). 

The complex pattern of calcium concentration can be measured by grouping or counting the 
number of cells for a given value of concentration. The results are plotted in Figure 3b. It is clearly 
seen that there exists a power law in the distribution of the number of cells (n) versus discrete calcium 
concentration (u). A least-square fitting of n oc u~'^, leads to the exponent of 7 = 1.26 ± 0.02. This 
implies the self-organized criticality in the complex calcium patterns. 

This result may have important implications to the calcium transport mechanism. Although 
the detail of intracellular calcium oscillation and communication is not clear, it is likely that the 
intracellular and intracellular interactions mainly occur locally. Thus reaction-diffusion dominates 
the process without much contribution from the convection mechanism. This is consistent with the 
physiological aspects of calcium transport and functioning [4,10]. 

5 Conclusions 

The finite state stochastic cellular automata have been formulated to simulate the reaction-diffusion 
systems of nonlinear calcium oscillations and waves using the transition rules derived from the partial 



differential equations. By using the proper stochastic and transition rules between different states, 
the finite state automaton can simulate the complexity of calcium transport. Numerical experiments 
show that complex patterns can arise from the initial random configuration due to the local transition 
rules between entities of certain nearest neighborhood and the details of initial configuration is not 
important. The power- law relationship between number of cells and calcium concentrations implies 
self-organized criticality in the complex patterns. 
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